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PREFACE 


Two feedback control schemes which maximize a terminal quantity while 
satisfying specified terminal conditions are discussed and compared. The 
schemes зге based on a linear perturbation from a nominal non-optimal path 
which does not, in general, satisfy the terminal conditions: The E 
have been programmed in the ALGOL computer language for evaluation and the 
programs are included in the appendices. 

I wovid like to express my appreciation to my advisor, Dr. Benjamin 
O. Lange, Dr. Richard Rosenbaum of the Lockheed Missile and Space Company, 
and Mr. Thomas Bullock, for their guidance and assistance during the 
preparation of this thesis. 
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I. INTRODUCTION 


In recent years several authors have treated the problem of deter- 
mining optimum control programs for nonlinear systems with terminal 
eanstraints. These problems arise in the design of control systems and 
development of guidance laws where it is desired to determine, out of all 
possible time histories of the control variables, the one control history 
that maximizes (or minimizes) one terminal quantity or cost function while 
simultaneously yielaing specified values of certain other terminal 
quantities. 

The steepest-ascent or gradient method developed by eee Bry on 
and Nao which is a systematic and rapid numerical procedure, has 
proved to be successful in solving this class of problems. Improvement 
in the convergence time of the iterative process involved has been achieved 
by E by a method based on the earlier work of Bryson and DS 

Another successful method developed by Breakwell, ae and modified 
by Bullock® is a second variation method in the Calculus of Variations. 

The principal objectives of this thesis are to develop a simpler 
steepest-ascent program which will be understandable to the control 
engineer without a background in the Calculus of Variations and to compare 
the results and speed of convergence with the method developed by Bullock. 
In this way it is expected that the steepest-ascent program will prove to 
be a useful instrument in education and research, while at the same time 
through the comparison, illustrate the advantages and disadvantages of the 
two approaches to the problem. 

The method used in developing the steepest-ascent program is 


essentially a variation of the methods of Bryson and poniendo. and 
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Rosenbaum , hence it is restricted to problems in wnich only the deviations 
in the control variables and adjustable parameters are considered in the 
performance index. It is also restricted to problems in which the pay-off 
quantity is a function of the terminal value of the states. This varia- 
tion includes a terminal-error control scheme which maintains a bound on 
the terminal constraint errors, hence reducing the total number of itera- 
tions required to converge to the optimum since larger deviations from the 
nominal trajectory can be tolerated while still meeting the desired terminal 
conditions: 

A numerical example is given of a rocket ascent trajectory into a 
circular orbit of maximum altitude. Provision is made for a two-stage 


rocket with optimization of the inter-stage coast duration. 





Il. STATEMENT OF THE PROBLEM 


Given a system which can be described by a set of non-linear (or 
linear) first order, ordinary differential equations, determine a control 


=t = 7, to maximize 


history u(t), in the interval to 


p = JM (1) 


subject to constraints 


ЕИ 0, (2) 
Ра : 
— e fé[x(t),u(t),t], (5) 
ОВЕ 
tg; T, and x(to) given; (4) 
where 
v (+) 
° an mx 1 matrix of 
ut) = j ” control variables, (5) 
a(t) 
x, (t) . 
` an a marco cuate 
x(t) = А аттар ее; (6) 
ЖО 
in а а х 1 matrix of terminal 
y | constraint functions, each (7) 
5 | ^ of which is a known function 
те on ZUR 
а 
f : 
118 ап п х 1] matrix of known 
| К | AM ult), (8) 
u j le ss ten 
2 equations; 
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Ø is the pay-off quantity and is one of the 
states, namely x, (T); (9) 


T is the terminal value of the independent variable. 


ОМ асп iteration iL с cesi vec. го minimize tie mean value of a 


positive definite quadratic form in the control variable deviations: 


с = Ih Se) une (10) 


where the superscript (') indicates the transpose of a matrix, and би(т) 
are small deviations in the control history from a nominal non-optimum 


Urscectory. 





111. STEEPEST-ASCENT METHOD OF SOLUTION 


A. BACKGROUND 


Bryson and Denham, in Ref. Z, considered terminal control of non-linear 
(and linear) systems for minimum mean values of a positive definite quad- 
ratic form in the control variable deviations. That is, it was assumed 
that a nominal control history had been determined which caused the vehicle 
to arrive at the terminal point with desired values of certain specified 
terminal conditions. Small deviations from this nominal trajectory were 
considered which might be caused by disturbances, inaccuracies in the · 
Gata, inaccuracies in the control system, etc. The problem was to determine 
small deviations from the nominal control so that the terminal constraints 
would be satisfied in spite of the disturbances. 

In the present paper the nominal trajectory is determined by guessing 
a reasonable control variable Dc For example, in & rocket trajectory 
problem one might choose an initial launch angle and & gravity turn with 
zero tnrust angle throughout as is done in the numerical example. > 
more, it is desired not only to determine control deviations which result 
in meeting the terminal constraints, but also to maximize the terminal 
value of one of the states while minimizing & performance index. This 
optimization scheme is a variation of the so-called Lambda Matrix Control 


4 3 
feedback method described in Ref. Z and the convergence method of Ref. 2. 


B. DERIVATION OF EQUATIONS 


The optimum programming problem can be solved systematicaily and 
rapidly on a high speed digital computer using the steepest-ascent 
technique. As stated in (A), this technique .starts by guessing a nominal 


control variable program, u*(t), and solving the set of differential 
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equations (3) with initial conditions (4), to determine a nominal tra- 
jectory. This trajectory, in general, will neither maximize ф nor will it 
satisfy the terminal constraints (2). 

Consider small perturbations in the control variables, u(t), about 


the nominal, where 
E) u) u*(t) (221) 
(The superscript (*) indicates terms evaluated along the nominal trajectory.) 


These perturbations will cause perturbations in the state variables, öx(t), 


where 


et) = x(t) - x*(t) (12) 


Substituting these relations into the system differential equations (3) and 


expanding f in & Taylor series about the nominal the result is, to first 


order 
<~ (8x) = F(t)öx(t) + D(t)su(t) (13) 
where 
an n xin matrix or 
R(t) = E partial derivatives 
— Vox/ ? evaluated along the 
nominal trajectory; 
and 
all HB x Mm matrik of 
p(t) = E partial derivatives 
"^ iQu/ ? evaluated along the 


nominal trajectory 


To determine the effects upon the terminal conditions f and y we 


introduce the linear differential equations adjoint to (12) defined as 


аф 
ар (7,5) - -Ф(Т,5)ғ(%) (14) 





where 0 is an nx n fundamental or state transistion matrix whose elements 
give the sensitivities of the terminal states to perturbations, 5x(t), 
alene the trajectory. (See Ref. 7). Initial conditions for these equations 


are specified at the terminal time, i.e., 
o(T,T) = I, the identity matrix (15) 


hence numerical integration proceeds backward in time. 
The solution to (14) provides a solution to the linear perturbation 


equations (12) at the terminal point: 


a) ax (to) qe 9(T,7)D(r)5u(r)ar (16) 
t 
Since f and y depend on the terminal values of the states, small 
deviations, df and dy may be calculated from the solution to (16). Con- 
versely, if 5x(T) is known by specifying values of 40 and dy, the corres- 
ponding control history may be calculated, which is what is done. Re- 
writiwe (L6) 


ex(T) — e(T,t)ex(t) fe o(T,r)D(r)8u(lr)är = O (17) 
t 


In order to minimize the performance index subject to the constraints 
(16), the method of Lagrange multipliers is employed. Multiplying (16) by 


a matrix of Lagrange multipliers 


pee 5 ; (18) 





(where (16) is written to include only those states which appear іп 2 ала 


y, hence q-1l equations) and adjoining the result to (10) 


ТЕ 
MEU ОСЫ) Jd [5 Su'(r)8u(7) -v'0(2,7)D(7)éu()] dt 
(19) 


Completing the square 
Tri 
E I I) oO) ax (e) ] S E (Su(r)D'(1)0'(T,7)v) '(Su(r)-D '(7)9 '(T,7)v) - 
t 
- 3 v'0(1,7)D(7)D'(7)o*(T,7)v |ar (20) 
To minimize C subject to the control variations choose 
зи(т) = D'(T)o '(T,T)v (21) 


Substituting this relation in (17) in order to find v 


р 
sums top) ex(t) -/ Ф(Т,т)р(т)0'(т)Ф(Т,т)удт = O (22) 
t 
Define the controllability matrix 


Ј = uf aaa me tener (23) 


Note: This integration may be performed simultaneously with (14) by 


numerically solving 


= = -S(T,t)D(t)D'(t)o '(T,t) (24) 


ILD а а сова ново 
ТӨШ) - 6 (25) 


Equation (22) may be written 





ои о јеж о) = ЈУ (26) 
solange LOr V 
yv = J [8x(t) - Ф(т,6)5х(%)) (27) 


Substituting this relation in (21) 


Su(t) = D'(t)o'(T,t)J "[sx(T) - 6(T,t)8x(t)] (28) 


which is the perturbation in the control history which satisfies the con- 
straint (17) while minimizing the performance index (10). 
As stated earlier, 5x(T) is determined from the specified values of 


dp and dy where 


do ex, (T) (29) 


and 


ay = vlex(T)] | (30) 


As yet, nothing has been said as to how one chooses the desired pay- 
off improvement, dq, or the desired improvement in meeting the terminal 


constraints, dy. The latter is normally chosen as 
ay = -y | (51) 


that is, the negative of the total error on any iteration is chosen as 
the desired correction specified on the following iteration. The problem 
of specifying dp is more complex and is the subject of the next section. 
It is worthy to note at this point that the Lagrange multipliers, 
which are error feedback terms, need not be computed at every point on 
the 2. Sufficient accuracy can be obtained in computing the 
control deviations (21) by solving (27) at discrete intervals and using 


the result until the next "sampling time". This reduces the number of 
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times the controllability matrix must be inverted per iteration and materi- 
ally improves the running time of the program.  Experimentation will re- 
veal how large the sampling interval can be made. Since the controllability 
matrix is singular at the terminal time, T, new values of the Lagrange 


multipliers should not be computed too close to the end of the trajectory. 


C. METHOD OF SPECIFYING THE IMPROVEMENT IN PAY-OFF 


In general, one does not know how far from the optimum a given nominal 
trajectory will be. It is, therefore, difficult to guess how much pay-off 
improvement to specify initially. However, it is possible to compute & 
value of ap which will result in a trajectory that satisfies the terminal 
constraints. This is done as follows: 

The changes in the control variables required to meet the terminal 


constraints with the pay-off unconstrained are given by 


sult) elf. 


al 
A zd dy (52) 
0 


where 


A = O0(T,t)D(t), without row 1 or colum 1. 


This is the same as the basic control equation (28) with ex (to) equal to 
zero and &x(T) containing only the terminal constraint terms. 
The change in pay-off, dd, that will be produced by a given change 


in the control variables is, from (16) 


do = ex, (T) Jf А ба(т)дт (52) 
0 


where A, is the first row of the A matrix. Substituting (se) in (32) 
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where A, is the first row of the A matrix. Substituting ет (52) 


an 5x. (T) = E Az "e J A vw dy (34) 


0 to 


Fouation (34) gives the change in pay-off associated with adjusting 
the control in order to meet nel constraints. "This value is used 
on the first iteration. 

Equation (3h) is also used to compute a value of the pay-off cor- 
rected for terminal errors, i.e., Ọ + dọ is the value the pay-off would 
have achieved had the terminal errors been zero. Hence, one can determine 
whether an improvement in the pay-off was ОО achieved or if ап 
apparent improvement was a result of larger terminal constraint errors. 

On subsequent iterations, one of three methods is used to compute do. 
A value equal to 25 per cent of the nominal value of E computed and 
stored. This quantity is called dq** and is used in method (2) below. It 
is a fairly arbitrery choice but should be made as large as seems reason- 
able. The program will automatically adjust it if it is too large. 

Method 1. Choose dO to satisfy the terminal constraints 

with the pay-off unconstrained as described above. 


Method 2. If | av | S €, where & is chosen as reasonable tolerance on 


the terminal constraints then 





a (35) 


21. 

where i is a count of the number of times method 
(3) has failed. This has the effect of halving 
the improvement specified each time a run is 
unsuccessful in improving the terminal errors or 


the corrected value of the pay-off. The program 


LL 





terminates when, while executing this method, 
dp becomes less than a pre-set number. A 


final run is then made using method (1). 


Mechod ou lT Lay) > e the following questions are asked: 
Were the errors on the current iteration 
smaller than on the preceeding iteration? 

Was there an improvement in the corrected value 
of the pay-off? If the answer to either 

que tilon is no, the mom as considered unsuc- 
cessful, the control history is replaced by the 
previous control history, and method (1) is 
used. If the answer to either question is yes, 
dp is set equal to zero and an attempt is made 
to satisfy || < <. = this test fails a 


second time, method (1) is used. 


D. TWO-STAGE ROCKET TRAJECTORY WITH COAST PERIOD OPTIMIZATION 


In many Orbit injection applications, such as the Gemini-Titan LI 
system, the launch vehicle is made up of two powere stages. It is there- 
fore of interest to consider the effect of an interstage coast phase on 
the maximum altitude obtainable. In this section a method of calculating 
the optimum coast duration is derived. 

The basic equations, (1) through (12), are the same. The linear 


perturbation equations (13) may be written 


а [5х(+)] = F(t)8(t) - D(t)8u(t) 4 B(t)8c (36) 


where 


antn xema rix Of partial 
B(t) = (3% derivatives of f with 
5% оре он ве сена о 
duration, c. | 


(37) 


T2 





The sclution to (37) is 


Sam) - 2/5 /)5х(5) SEL CDU OE + `®(т,т)в(т)всат (55) 
i 


t 


Tne performance index becomes 
lide 1 2 | 
E .if би'(т)би(т)ат + zd Se l (39) 
16 


where A is simply a weighting factor. 

Before attempting to minimize (39) subject to the constraint (38), 
a method must be derived to evaluate the last tem in (38) which is the 
change in the terminal values of the states due to a change in the coast 
duration. Since с is an adjustable parameter which dces not appear ex- 
plicitly in f, the partial derivatives cannot be evaluated directly. 


However the desired term may be calculated by the following method: 


Define 
Ui = Stage 1 burnout time 
5 = Stage II ignition time 
hence 
е = ОИ (40) 


Since бс is a small time increment we may write 


to+At 


- xft = f dx 
x(t, + At) х(5,) i um 
2 


x(to)+Ax 


il 


x(to) 


| 


(x(t,) + Ax) - x(t5); 
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consequently we may make the approximation 
LONE CU 
x(t, * 8c) - x(t,) 3 (2), 5c (41) 


Now define x, as the states evaluated with the thrust off (uncontrolled), 


and x, as the states evaluated with the thrust on (controlled). From (41) 


x(t, + ёс) x(t,) + X, (t,,)80 


= + Y 
x (t, + c) x(t,) ž (t,)50 
where (*) indicates differentiation with respect to time. Subtracting the 
above expressions we have 
x (to + dc) - x (t, + Sc) = (x, (t,) - X (t5)]6c 

Define 

cx, = [x (t5) - X(t.) Ide Ä (42) 
The quantity ox. is à perturbation in the states occurring at time 

= + 

E to OC 

Teee UP ori the thrust for a period 6c. The question remains: 


How does this perturbation propagate to the end of the trajectory? The 


answer is clearly 
&x(T) = 2(Т,%, + Sc)&x, (43) 


Finally, (38) may be written 
T 
=> o (Tux (a) +f Ф(Т,т)В(т)5а(т)ат з 2(Т,%, + Sc)öx_ (44) 
E 


where the last term is zero prior to time t5 NO 


1h 





Introducing the Lagrange multipliers, v, and adjoining (44) to (39) 


T 
(= E би'(т)би(т)ат + = Ме“ + v'Bx(T) - 5'o(T,t)J5x(t) (45) 
А | 


= 


T 
- vif (T,r)D(r)du(r)är - DEDUT t de Sc)Ex, 
ib 


As before, completing the squzre yields the optimum control change 
бщ(т) = D'(r)o'(T,T)v (46) 
Puga степ ла отп Lhe Optimum coast change is 


be = Fv! Olt, t, + бе){% (%) - % (6, [ (47) 
Substituting (46) and (47) into (44) 


T 
&x(T) = 9(Т,%)5х(%) Jd O(T,7r)D(r)D' (1) O'(T, 7) atv 


+ + у А T Se) [x (t,) - x (t,)1I[x, (t,) - ac QUARTO +5c)v 


А 
Define 
Douce os oT. text) (49) 
T 8 
Ј "i ee (50) 
ЕЕС 1 625) 
G= O(T,t, + 8c)A A' OT, t, + бе) (52) 


Rewriting (48) 


А = Ју + = бу (55) 
eal | : 
у= (2+8) A (54) 


15 





Finally, substituting in (46) and (47), the control and coast variations are 


5u(t) - D'(v)e'(7,9) (3 i M" ^ (55) 
56 - SE A'O(T,t, +80) (9 «з [öx(T) - 9(T,t, +5c)5x(t, +50)] (56) 


Since the last term in (44) is zero prior to time t, + 5,2 the term 


2 
G/A in (55) and (56) is also zero prior to that time. 


Equation (56) should be evaluated at t = t, - 8c, but since 8c is the 


e 


unknown, it is evalvated at t This does not introduce an appreciable 


5° 
emrer if бо às small. 

In some numerical integration procedures the terminal value of the 
independent variable cannot be changed once the integration has begun, 
hence the change in coast time cannot be added immediately. This problem 
is solved by evaluating (56) on each forward integration (except the nominal) 


а OC Ё O, re-integrating the latter portion of the trajectory from 


54 Tow. 


E. COMPUTATIONAL PROCEDURES 


As stated earlier, this steepest-ascent technique requires the use of 
a high speed digital computer. The sequence of operations is summarized 
here. 

l. Compute the nominal path by integrating the system differential 
equations (5) with & nominal control history and appropriate initial con- 
ditions and store the time history of the state variables at reasonably 
smell intervals. Print out the values of © and ү. 

2. Integrate the adjoint differential equations (14) backward, 
evaluating the partial derivatives on the nominal path by reference to 
the states stored in step (1). Simultaneously integrate the controllability 


matrix equations (24). Store the results at the same interval as the states. 
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3. Select desired terminal condition changes, dọ and dy, as ex- 
plained in Sections B and C. 

4. Compute and use the new control history while integrating the 
system differential equations forward. Again print out the values of Фф 
and v, unless the next step applies. 

5. If the two-stage rocket problem is being solved, compute the new 
coast period in step (4). Transfer the 5200270 of the second 
stage control history to correspond with the new coast time. Integrate the 
system equations from 54 to the new terminal time. Print out the values of 
Ф апа ү. 

6. Repeat procedures (2) through (5) until the pay-off improvement 
in step (3) is less than a preset value. At this point, use method (1) 
described in Section C to select dp and complete step (4) and (5). This 
has the effect of eliminating any remaining errors in the terminal 
constraints. 

Т. Punch cards or store the control history on tape and terminate 
the program. 

Before concluding this section, a few general factors of great 
importance in this type of numerical calculation should be discussed. 
| Ihe programmer must exercise great care when working with values of 
type real (or floating point). Often a calculation is made where the 
result is expected to be an integral value such as 4/2 = 2.000 ..., 
however, due to the binary, octal, and decimal conversions which take 
place within the computer, а may come out 1.9999...99. This 
problem occurs when trying to generate array storage indices based on a 


value of the independent variable which is a floating point quantity. 
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IV. SECOND VARTATION METHOD OF SOLUTION 
A. OUTLINE OF THE METHOD 


Bullock? has derived a feedback control scheme based on the second 
order variational theory in the Calculus of Variations. The method is 
outlineä here in sufficient detail to solve the problem stated in II. 


The differential equations to be satisfied are 


ш t) (57)‏ )ا 
A (EA (58)‏ 
M F -Q M‏ 
o»‏ )50 4:0 
N ON‏ 
and‏ 
b = M'a - N'e | (60)‏ 


in order to maximize 

puso. T) (61) 
and satisfy the terminal constraints 

ү = y[x(T), T] 9 
Define the variational Hamiltonian 


H = A'f (65) 


The elements of Eqs. (59) and (60) are 


E Oia | (6%) 
X OOA 
= : 
Жш | (65) 
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S - Н -Н H H' (66) 
Der ЖАЛ Дй 
= i 
en (67) 
= | 
a= He HH (68) 
xu uù u 


where the subscripts indicate partial differentiation in the usual sense. 
The initial conditions for (57) are given and the terminal conditions 
262 (58), (59), and (60) are 


AM - vw) (69) 
5-І 


where the components of the column vector v are sensitivities of the pay- 


off Y to changes in the terminal constraints y; 


wMT)-ri-wQG v) oc (70) 


RR, X 
where I is the identity matrix; 


N(T) = -yty (71) 


MU 


H 


b(T) = -y, dy (72) 
where dy = -y. 


The perturbations in the control history are given by 


EM C ae _ _wm-l t 
Su(t) = u u = Ka + Hs + ft SA) (73) 


where 


ЗА = (MY) "(max + by. (74) 


In order to minimize a performance index 
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Q 
i 
mir 


T 
J S6u'(r)8u(7r)aàT oy 


E in the above equations is modified by adding an arbitrarily large 
negative constant, K, which is reduced in magnitude as the program con- 
verges. This has the effect of constraining the magnitude of би. 

Since the terminal conditions on the adjoint equations, A(t), depend 
initially on the choice of v, a method is given which will improve the 
accuracy of these terminal conditions on subsequent iterations. 

Equation (74) is an expression for 5A at any time, t, but since M is 
singular at T, it cannot be evaluated directly. However, if a point (t,) 
Ш- сповеп süfficiently far from this singularity, the following equation 


can be integrated from 54 cowl: 
A Ho (76) 


men 1nbial eenditiens 


3Mt,) « Qe) For öx + QM (77) 


The solution to (76) is then added to the current values of A(T) prior to 


the next backward integration of (58). 


B. COMPUTATIONAL PROCEDURES 


As in the steepest-ascent method, this method requires the use of a 
digital computer. The sequence of operations is summarized here. 
l. Select a nominal ТН history and initial values of v. This 
can be done by starting with a control program and v generated by the 


Steepest-ascent method, where 


2l 





mom. Ea. (55). 

2. Integrate the system differential equations as in the steepest- 
ascent method. 

3. Integrate Fas. (58), (59), and (60), backward with the appropriate 
terminal conditions. If the determinant of M changes sign or Hy becomes 
positive, store the current value of the time and stop the a. 

The reason for this is explained below. 

1. From Eq. (73), compute the new control history while integrating 
the system equations forward. 

>. Compare the value of ® and y obtained to those obtained on the 
nominal trajectory. If the pay-off, 9, or the terminal constraints, y, 
have become worse the run is considered unsuccessful and a tighter bound is 
placed on ôu by increasing the magnitude of K. If, on the other hand, Ọ 
and y are the same or have improved, the run was successful and (3) and (4) 
are repeated. 

6. The program is terminated when no change occurs in the pay-off 
er the constraints and [x] << ІН 1. А this pointithe control history a 
stored on punched cards or tape. 

7. If in step (3) the determinant of M changed sign (this condition 
is called a "conjugate point" in the Calculus of Variations), or il. 
became positive (which MEL. the Legengre condition is not satisfied), 
the integration in (4) is begun at a slightly later time than this con- 
dition occurred. Normally, on subsequent iterations, this point will 
move backward to the beginning of the trajectory and disappear. 
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V. PROGRAM EXAMPLES 


A single-stage rocket trajectory problem as described below was 
programmed utilizing each of the methods discussed. The ALGOL computer 
language was used and the programs were run on a Burroughs B-5500 digital 
computer. 

Assuming the rocket is launched from an airless, non-rotating Earth, 


the state equations are 


X, = Y = V sin (y) 


O - – cos (Y) 


W T 


A — 


Wo Isp 


E c Ds cos(u) Eins 
О О 1 t „2 


қ y a sin(u) |" cos U соз (y) 
2 Y 
T 


(V Wo T E 


pS 
= 
|| 

= 


Wo Isp 


where r = altitude measured from the center of the Earth, V = velocity, 

Ү = flight path angle, go = gravitational acceleration at the Earth's 
surface, T = thrust (assumed constant), u = thrust angle (measured from 

the velocity vector), Wo = initial weight, Isp = specific impulse, + = time, 
Н = universal gravitational constant, © = downrange angle. The initial 


eonditions are 


(= R, (Earth radius) 


.e(0) 2 0 
Ү(0) = 100 ft/sec 
v(0) = 89.87 degrees 
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It is desired to place the payload in a circular orbit of maximum altitude, 


hence 


= 
|| 


> = r(T) 


% t=T 


Appendices (A) and (B) contain listings of the steepest-ascent program 
and second variation program respectively. Comments are inserted at strategic 
points which explain the sequence of operations. 
The steepest-ascent program contains logic for a single or dual stage 
rocket. It was run in the single stage mode to generate a nominal control 
history for input to the second variation program ana to compare results. 
It was also run in the two-stage mode to test the coast optimization (овом 
The input data for the steepest-ascent program are 
l. Initial velocity (feet/secona) (must be non-zero). 
2. Launch angle (degrees). 
3. Duration of the first stage burn (seconds), for single stage 
Meckets this quantity 15 the Total burn time. 

|. First stage thrust (pounds). 

2. Second stage thrust (pounds), for single stage rockets zero 
Ts Input. 


6. First stage fuel flow rate (pounds/seconä). 


1. Second stage fuel flow rate (pounds/second), for single stage, 
any non-zero number. 


8. Rocket liftoff weight (pounds). 
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ШЕ 


Tis 


125 


Hc 


Second stage weight after separation (pounds), for single 
stage, any non-zero number. 

Initial value of coast duration (seconds), for single stage, 
Zero. 

Duration of the second stage burn (seconds), for single stage, 
zero. 

Coast weighting factor, A. 


Number of stages (1 or 2). 


The input data for the second variation program are: 


T 


2. 


9. 


Initial velocity (feet/second). 

Launch angle (degrees). 

Duration of rocket burn (seconds). 

Ihrust divided by initial weight (pound/pound). 
V 


ЈЕ 


vo 
Integration step size and data storage interval. 
K (See Section IV-A). 


kKormmnal control history- 


Other parameters which the user may desire to change must be changed in- 


side the program or incorporated into the READ statement- 


For the single-stage runs, the following input values were used: 


Initial velocity... «.o.o..o.o.. OO sec 
EE lenene ees O O E E rees 
. Duration of burn..... 2222. ‚220 seconds 
TEC teu. ОО е ЛЫ 0 БОРО 5 
а ООРУН 1433.3 pounds/ second 





Liftoff weight..... 


еееееее 
oe 9. е е ее ееге әееее 
еееее 


For the two-stage runs, data 
Initial velocity... 
launch &ngMe...... 
First stage burn... 
First stage thrust. 


Second stage thrust 


First stage fuel flow 


Second stage fuel flow. 


MCOT verent. eas. 


Second stage weight 


Cosa stsoauratıion 


Second stage burn.. 


ОАЕ ИЗБОРЕ ОТ ое 


Number of stages... 


cc Meu 0 Pounds 


2222222 


for the Gemini-Titan II system were used: 


.100 feet/second 


2..09.9( %о 98.05 деггеев 


see. aD) seconds 


(997588000 pounds 


. . . „ „100 ,000 роцла 


1666.6 pounds/ second 


DE] pounds/ second 


SL OOO БОШО 


../0,000 pounds 


Seconds 


.....180 seconds 


0.1 


2 
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VI. RESULTS AND DISCUSSION 


A. STEEPEST-ASCENT VS. SECOND VARIATION 


When the steepest-ascent program was run in the single-stage mode, 
the nominal en, Mm eure or ТОСО Е The errors 
in meeting the terminal constraints on the velocity and flight path angle 
were 426 feet per second and 1.156 degrees. On the third iteration the 
altitude was improved to 260,427 feet with terminal errors of 0.67 feet 
per second and 0.025 degrees. The program was terminated when the desired 
altitude change, аФ, became less than 5,000 feet. At this point fifteen 
iterations had been completed. The terminal altitude achieved was 318,126 
feet with terminal errors of 0.64 feet per second and 0.003 degrees. The 
associated control history was punched on cards and values of Ya and Vo 
were printed out. The program ran Ter minutes and nine seconds. It is 
estimated that this time would be about halved if the program were run on 
an IBM 7090 computer. 

The output generated in the steepest-ascent program was used as input 
to the second variation program. As expected, the nominal trajectory 
attained an altitude 318,126 feet. On the succeeding forward integration 
the trajectory was totally unreasonable. The control deviations were 
made smaller by increasing the initial magnitude of K but this failed to 
improve the results. Small variations were made in the input values of. v 
which caused relatively large changes in the results but it was not clear 
how to make adjustments which would improve the performance of the program. 
The eh time was far in excess of the steepest-ascent program, taking 


over three minutes to compile and complete just one iteration. 
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B. TWO STAGE ROCKET TRAJECTORY WITH COAST OPTIMIZATION 


As stated in Section IV, the input data for this problem were those of 
the Gemini-Titan II launch system with an arbitrary choice of ten seconds 
for the initial coast duration. Lambda, the coast weighting factor, was 
chosen as 0.1 since earlier runs indicated that a value of 1.0 caused the 
coast variations to be insignificant. This choice proved to be satis- 
factory. 

The nominal trajectory attained an altitude of 392,564 feet with 
terminal errors of 25 feet per second in velocity and 0.84 degrees in flight 
path angle. On the first iteration, where the program attempts only to 
meet terminal constraints, the altitude was improved by 18,018 feet, the 
terminal errors were 0.09 feet per second and 0.0025 degrees. On the 
fourth iteration the coast time was reduced to eight seconds. At this 
point, a new nominal for the portion of the trajectory following first 
stage burnout was computed using the new coast period. The result of 
this change in coast was an improvement in the velocity constraint of 
12 feet per second, and a degradation of the flight path angle by 0.25 degrees. 
The terminal altitude achieved on this iteration was 573,450 feet. On the 
fifth iteration the coast period was reduced to two seconds, this accounted 
for an altitude improvement of 30,973 feet. On the sixth iteration the 
coast period was reduced to zero, this improved the terminal altitude by 
13,477 feet. In both instances cited above where the coast period was 
changed, the terminal constraint errors were diminished. 

The altitude improvement specified on the second iteration was 50,000 
feet. The program was allowed to run until this figure was reduced to less 


than 1000 feet. This proved to be rather wasteful as the terminal altitude 
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failed to improve significantly after the desired improvement was reduced 
to 6,250 feet, which occurred on the eleventh TM In Fig. 1 the 
terminal altitude, corrected for terminal constraint errors, achieved on 
each iteration is illustrated. It clearly shows that twelve iterations 
were sufficient to converge to the optimun. 

Figure 2 illustrates the convergence of the control history. The 
discontinuities which occurred are due to using discrete feedback at 
twenty second intervals rather than continuous feedback. The curves are 
of different lengths due to the changes in coast period which occurred. 

Figure 3 illustrates the action of the terminal error control scheme. 
Each time the errors became excessive they were reduced to essentially 
zero in one iteration. 

The last terminal altitude achieved was 616,573, a 57.2 per-cent 
improvement over the initial nominal. 


C. SENSITIVITY OF THE STEEPEST-ASCENT PROGRAM TO CHANGES IN INITIAL 
CONDITIONS 


The sensitivity of the steepest-ascent program was tested to determine 
the capabilities of the convergence scheme. Lauch angles of 89.07 and 89.95 
degrees were tested on the single stage trajectory. The resulting terminal 
altitudes were 196,015 and 803,428 feet. In the latter case the terminal 
errors were 1858 feet per second in velocity, and 40.4 degrees in flight 
path angle. In both cases the program converged in twelve iterations to 
approximately the same optimum altitude. 

An attempt was made to Ets the problem given a ninety degree 
launch angle which failed. Further testing at lower launch angles was not 


eccomplisned due to computer time limitations, however, it is believed that 
the tests conducted amply illustrate the virtue of the method. 
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VII. CONCLUSIONS 


n 


The terminal control feedback scheme due to Bryson and Benham and 
the method due to — leaving the pay-off unconstrained in order 
to satisfy terminal constraints have been combined and altered as necessary 
to produce a simplified steepest-ascent optimization procedure. This 
procedure has been shown to be successful in solving a typical rocket 
trajectory problem including optimization of an interstage coast period. 

In this closed-loop procedure the change in control is computed at 
each point utilizing continuous or discrete information on the deviation 
from the previous trajectory. It is closed loop because the program 
continuously (or periodically) checks on how it is doing in its attempt 
to satisfy the terminal constraints. The advantage of this procedure is 
that larger deviations from the nominal can be tolerated while still 
maintaining a bound on the terminal errors. It is, therefore, possible 
to move rapidly toward the optimum trajectory as is illustrated in Fig. 2. 

The second variation method due to Bullock? has been applied to the 
same problem with questionable results. Consultation between the author 
and Mr. Bullock has failed to uncover possible flaws in the theory or the 
programming technique. Bullock has shown in several examples that the 
method is successful, however, in each example the Hamiltonian did not 
depend explicitly on time. In the rocket problem this is not the case. 
Although this fact has no theoretical bearing on the problem, it is the 
only major difference between Bullock's examples and this problem. 

Experimentation revealed that the rocket problem is extremely sen- 
sitive to the choice of v. This was not the case in the examples 


presented by Bullock. It was thought that the values of V generated by 
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the steepest-ascent program would be very close to correct, however, in 
view of the failure, this premise was laid open to question. Hence, one 
reasonable conclusion that can be drawn is that v was chosen incorrectly 
and that there is no intuitive or analytical method presently available 
to make the proper choice. There exists, of course, the possibility of 
error, but considerable time and painstaking effort has been expended to 
minimize this possibility. 

Failure of the second-variation method notwithstanding, some con- 
clusions may be drawn with respect to the advantages and disadvantages 
of the two methods. 

l. Understanding the theory involved in the second-variation method 
requires a background in the Calculus of Variations, whereas the steepest- 
ascent method presented does not. 

2. The second variation method apparently requires a good estimate 
of the sensitivities, V, while the steepest-ascent method only requires a 
guess of the nominal control and will tolerate a fairly poor guess. 

5. The second variation method requires backward integration of the 
2n? * 2n equations M, N, b, and A (where n is the number of state variables) 
and a matrix inversion at every step of the forward integration. The 
steepest-ascent method requires backward integration of 6 and J which is 
less than 2nf Bene since J is symmetric. The matrix inversion can 
be done at less frequent intervals using the sampled data feedback method 
suggested. The second variation method thus requires more computer time 


per iteration. 
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VIII. ADAPTATION OF THE PROGRAMS TO OTHER TYPE PROBLEMS 


An attempt was made to generalize each program so that they could be 
easily adapted to other problems. This required a large number of sub- 
scripted variables and matrix multiplication loops in the subprograms 
containing the differential equations. Since the subprograms are called 
twice per integration step they must be as efficient as possible. The use 
of subscripted variables and loops is most inefficient and results in more 
than doubling the running time. 

Adaptation to another problem is still relatively simple however. 

The programs contain sufficient comments to indicate where necessary 


changes must be made. 
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APPENDIX A.  STEEPEST-ASCENT COMPUTER PROGRAM 





PEGIN COMMENT ERNEST Ce LUDERS ROX 219 STEEPEST-ASCENTS 


INTEGER FLAG^FAIL»CUUNT^QUITJ 

АЕА Т1», СОАѕЅТ У, РЕ РЕ 2» 15Р1», 15Р2» НЕЇ» МЕ М1, 25 Т2 ММ1, ММ2 
RE» MASnRAT, TUVERWO, GO, X» ISP, TF» HH» SAMPLETIME,» 
VÜ»PREDHFÜ'PREHFO»DHFC^HFO» J1 2» JI 34 PREVF »PREGAMF 
PREHF» VF» DGANF» 3Q0UND^, BCUND2» DET» EX» СМК» Р» 


MM3,LAMDAS SUM, ULDCOASTo TB2» Go DHF 5 


INTEGER MAXTHOEX»L»M»]> ITER»STAGES»OLDMAXINDEX» 


COMPUTECUASTsCUASTCUMPUTE 3 
Renee tle cts SJ» AC1:31])» YP£O0:4,0:2:110)5 


SAVE ARRAY XP,0LOXPEO:4,0:200)]5»LP(0:2250:200)» ERRIO? 3)» 
ЕРО ЧИО БА Са СПКГІ 15371; 
[VWP[1:2»1:21]» PH] 5,0LCPHI(O?:200)»JEI» 


ИИ DISS 
LABEL lis, L2» | 3» 4» 5» 6»нЕЦ 2 


ТС ПИНКИ РАНЕ БИ АДИ IA THE DETERMI", 
NANI UF J = j | 
ВИИИ /)/» X33» 3223.11» //» X18» "CHECK MATRIX s ", 
ОЕЕО 27333, SECS, 11s. 775 
F2 (X20 "THIS RUN REQUESTED AN ADDITIONAL: "s ИИ)» 
fe uM» " EFT OF UDE aT / 7/33 X30 F92» 
М meee, seo ОГ VELGCITY" o 7/75» X325» F11.6» 
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DE 5 ОР FLIGHT BATH ANGLES 4), 
A EN. ACPUALEY YIELDED AN AODITICNAL:;"s //)» 


omen Os TERMINAL ACTITUDE ="5 Fi0.2° 
" FEET", /» x20, "TERMINAL VELCCITY"^X175 "s", 
MONIO ин FEET/SEC"» /» X20, 
"TERMINAL FLIGHT PATH ANGLE ME GO 


% DEGREES", /9 


OS ITIL VELOCITIES ALTITUDE = ог X157:9.25 


4 EEE SEC"), 
Fé (X20»"COAST DUKATION",X2C, "2" ,F10,2," SECONDS"), 
Р? (X20»"CORRECTED IERMINAL ALTITUDE =")F10.29X7>» 
"РЕЕТ"), 


NTO ОЕМ 3 2 ) 5 
DEFINE LOOPL = FUR Lels253334 DO- 1; 
ШЕК COT 0 (2»10)» 


COMMENT MATRIX MULTIPLICATION PROCEDURE GOES HERES 


COMMENT  FATRIX INVERSION PROCEDURE GOES НЕВЕ? 


CU se) Sms PRO€CNUtc CONTALNS THE SYSTEM DIFFERENTIAL 
EQUATIONS AND IS USED UN THE FIRST NOMINAL» ON RUNS IN 
WHICH THE CUAST TIME HAS BEEN CHANGED » AND ON THE FINAL 


PRECISION RUN5 





БРО БАКЕ FUENCTCT^ X» F3 

VALUE T» REAL T5 REAL ARRAY X» Ft1)]; 

PEGIN ТЕЛЕН, V» S» C^EFEE»TS; INTEGER |» у; 
CUE O IDET ATE IN THE CONTROL ARRAY FOR THE 
PROPER VALUE; 

AE RCA J E TE IRA ENDE X THEN I*i ELSE Is 
FCO A A С ТАИ ТИРА E Y er 
EA EID 
S e SINC(OXU410» C e—QUSCX[(4105 V e X[3]} 
је | = TIME TIME FROM STAGE TWO IGNITION 
BT < OR STAGES ENG 
TüVvER"UücFF1/n15 
MASRATe1-TOVERWOXT/[SP1J 
END ELSE TUVERWOEO3> 
IF 1212 AND T2ZTF THEN BEGIN 
TÜüvER"nUecFF2/n25 
MASGATe1-TOVERWOXIS/ISP23 ENDJ 
COMMENT THESE ARE THE SYSTEM DIFFERENTIAL EQUATIONS; 
BI) TONS 
Ғ(21 e VxC/nj 
ЕЕЗ] e (Q € GOXTOVERWO/MASRATIXCOSCFEE) = KxXS/R*23 
ВАМИ ССПа + OXS EN CREETAVO + Piel; 


ERNEUERT, 


COMMENT THIS IS THE BACKWARD INTEGRATION PROCEDURE WHICH 
CONTAINS THE ADJOINT CEQUATIONS ANO THE CONTROLLABILITY 


MATRIX EQUATIONS; 
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BRBCEDURE UIA Ств LS» LF), VALUE TB5 

REAL [85 REAL ARRAY LS» LFt1); 

ON A UR SR CL I, Sfo 1, 12, і3, (4, 1МҺ12 
INTEGER TJS 
REAL TS; 
КЕА 15565,75 3 М1. М2, МЗ, МД у М5, Мб, М7, М8 3 

REAL Т5М,ТСМ»М1»Л2,М3»ҒЕҒ»(13>114>133>134»143>144); 

COMMENT JT IS BACKWARD RUNNING TIMES 

Гете, | seen enc TAH; 

LR Loomer І 76-05 

СИЕ ПРУ А МАЗАРИ ОВ РАНА У. NTE ПА А ЛИ 

ТЕС 1265221» L3eLS[3); ле 541 

Bc pn об о је Пеле о Сиј; Léelstads 

Һ16:5(9); 72126151101) М836/5(111) NA4«eLS(1215; 

ІЗЕЕӘГІЗ)», Мберзімллі»; DZELSt 1517 N8eLSC15)]; 
CU] INTERBOLATE. FÜR ТНЕ. РКСРЕВАМА ЏЕС ПР HAE STATES 

R e CASAR» 11) + INTXCXP(193] *-N) + RES 

үле Сасіргз,і1) * INTXCXP(3»J] ч 09; 

S € SINCcDeXxPLS,I)I + INTXCXPL4>J) = Q))3 

ccc cusccQexPDL4.s110 ғ ТҮТХСХРІ8,М1 ” 099; 

SE € о је РАДЕ + ANT CAS = QF); 
TSeT-T12> CUMMENT TIME TU GY UNTIL STAGE TWO IGNITIONS 
IF T272 AND T2é4TF THEN 4EGIN 

TOVER^OcFF2/u25 | 
MASRATe4-* TOVERWOXTS/ISP2J END 
FLSE TUVERWO€OS5 


E TU ON STAGES 51 THEN BEGIN 
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TOVERWO«FF1/515 
MASRATe€1-TOVERWOXT/ISP1; 
ЕМО; 
COMME THESE STE THE AUJOINT DIFFERENTIAL EQUATIONS; 
LFLI] eTVXCXL2/RA2 +2XKXSxXL3I/R*3 = CxX(V - 2xK/(VxXR))XL4/ 
ї*2;} | 
CEE Ec o 
LF[3] ¢ +SxXLi +CxXL2/R +Cx(1 + K/(V*2XR))xXL4/R = GOX 
TOVERWƏXSF/CV*ZXMASRAT)XL43 
A IS LA ALSO SINO KA CV ОУ 
xL4/n5 
A о A AS E E A ЕВ 
R*2j 
Epon 5 
CT TOOTS tCOXLO/R ECK EI a K7 CVX 2xR)IxLayk = GOx 
TÜVERNOXSF/CV*2ZxMASRAT2xL83 
UOTE = VXSXLO/R =кхохртукжкеа + SxX(V = K/(VXxR)) 
xl 3/33 
LF(9 ) € 7"VxCxN2/R*2 *2xKxSXN3/R«3 - Cx(V - 2XK/CVXRIIXNAZ 
HR*2j 
LECTO] emio; 
LFC[11] «© FSXN1 +CAN2/R РОАН У СРЕО О ЈАНИЋ = GOX 
TOVERHOXSF/( VZ XMASRAT)XNA3S 
ШЕПТЕ RS NRZ S x (VEAC XRD) 
«мдин 
БЕГІЗ! ӘБЕС /Кх2 ТКТ /ӘЮха3 т СХСУ - 2xK/CVXROOXNA8/ 
R*2j 


LF{14] e :057 


HO 





БЕГІС 25158 FCXNO/R ICx(]1 + K/CV*2xR))xNH8/R * GOX 
TOVERNOXSF/CV*2xMASRAT) XN; 


ГЕГ1Т6 UX CNS = VXSXN6S/N "KXCXN7/R«2 =~ Sx(V * K/CVXR)) 
x«N8/R; 
FEE © CQEPHICT IJ) + INTKXCPRICUJ-Q)3 
TSM * GOxTÜVERVO/MASRATXSIN(FEEDJ 
TCM TERME OAIMASRAT*XCOSCBEE/VJ 
MI PESAS Fr EIS S 
M2 € «М3ХТ5М * NüxT[ CH; 
M3 € -N7XTSH * N8XxTCH;j 
COVENT THE IRE THE CONTROLLABILITY MATRIX EQUATIONS; 
LELITT E MIL FB Te 8252: ІрГІ916 M3*2; 
LFr20Je MixM25 LFLI[21Je M2xM3* LF[22)Je М3ХМ1; 


END LINBAR» 


COMMENT THIS IS THE FORWARD INTEGRATION PROCEDURE WHICH 


COMPUTES IHE NEW CONTRUÜLS 


PROCEDUISESESINUPE Als FO VALUE Т) 

REAC ЕЕ АСКАТ ху РОХ 

BEGIN Әди, 5» C» TB» NDS BI INDOY EEC ITSM, TENNIS 

E.S HIS 157 | 
АНА нанио» [32 LA, TEOCLOSEEL5,L63 INTEGER I» Jj 
TSeT-1T2; CUMMENT TIME FRUM STAGE TWO IGNITIONS 
ПЕ <il O TAGES St LAEN GEGEN 
TÜVER4OÜ«FF1/u15 


MASAS LO NERNOAT/ ISP1) ISPCeTSP1j 
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END ELSE TOVERuOeO; 
IF T2712 AND T2X4TF THEN BEGIN 
TOVERNUEFF2/W2> ISPeISP23 
MASRATe41-TOVE 3HoOxTS3/ISP25; ENDS 
IF T1230UN0x.9992999999 UR COMPUTECNAST=1 THEN 


GEGIN 


PERIOD Is CUMPUTED> 


COMMENT THIS TS THE SAMPLED DATA FEEDBACK SECTION AND IS 
HIT AT ГНЕ SAMPLING INTERVAL AND AT T2 WHEN A NEW COAST 
ІР COMPUTECUASTA1 THEN 
BEGIN 
INDe€G80UND/HHs5 
dOüguNDega0UNDESAMPLETIME; 
veer, Se ree NT IER Cle ins SU UU COO To 
BNI NDEMAXINVEA*IND; SUM «03 
COMMENT SINCE J GETS SINGULAR NEAR TF» SKIP THIS SECTION 


CLOSE T IF; 


ІР BOUND 2 TF*30 THEN 90 TO TOOCLOSE; 

ARIS ОО ЕСТУ е ОСТУ Т ЗЕРЕ 5,58 IND Joe 
Cex. ООО OR COMPUTERS IT THER ОС ТУТИ 
LANDA ELSE 025 
JEILt>2JeJ1NVL1>21IeJEIL2>1JeJINVE2,LIeLPL20>BKINDJr+ 
2 | 2%, 759999 UR COMPUTESDAST=I ТНЕМ 0117274 | 
LAHWDA ELSE 925 

ЧОО ЈОЈ дз е see OUTAGE 5» о јереј» РКО ј+ 
mur Vee 799999 OR COMPUTEGOAST=!1 ТНЕК 0(2»31/ 
LAMDA ELSE 0)3 


eis ieJINVE ls 3J JEJEIL3 re LIEJSINV( Ss LIELP(22sBK INDI] + 
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CIF 17202392999 OR COMPUTECOASTST THENDE1>3)/ 
ANOS ESE O) 
ОСО ТЕС ША 5ТЕ 1 THEN GCOMPUTEGOAST€O; 
Ре ТЕСТ): 
1МУЕКТСЈІМУ, 3,1) 3 
[Р “ісі TAON JEGIN 
КЫТ СЕНЕ AMA TRA CES SINGULAR AT TIME = "s 
56223179 
COMMENT IF THE J MATRIX BECOMES SINGULAR TERMINATE THE RUN 
GU [0 HELL ENDS 
ODPL, NPIL] € ÜLDXPIL^INDIZ 
ШОРА suu vesc BA INDEX TL = pem» 
ИЕСІ 1 е ОЧЕ ниво Понте 02 
МООР SUN e SUN SC PUR ОРО ИУ Е ТЕ ИЕЛ); 
ЛЕРІ ée DYF * SUM; SUM € 05 
ПЕШИЕ О ОИЕ ОИЕ 12, ЕВКІМПІХ(ХІІІ * TEMPCLIOJ 
ОВ {3] € DGAME - SUM; SUM e 0; 
EAR e ERRE] GERRI oo, 
FUR Le1»25»3 JU FUR Me15,2»5»3 DO 
EKREL] © ERRCLI * JINV(L, HIXDELCMIS 
COMMENT ERR IS THE MATRIX OF LAGRANGE MULTIPLIERS» THIS 
IS THE END UF THE SAMPLED DATA FEEDBACK SECTIONS 
TUOCLOSE: END; 
COMMENT THE CONTRUL 15 CUMPUTED BETWEEN HERE AND L4» NÜUTE 
THAT RE ENDE SETSECH HP U PEDSONEESTEFFAHEAD OF THE CURRENT 
TIME TO PERMIT [NTERPOLATION; 
Li: IF T230UND2x,999999299 THEN 


BEGIN IF T0 THEN ВЕСІМШМей; GO TO 12 END; 
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(3: 1 е BNJUND2/HH + 15 —BOUND2 € BOUND2Z + НЫ; 
IM AS NEN GO PALAS 
2 НАИС РАТА е леза BKIND € MAXINDEX = T3; 
TS¢BOUND 23 
IF BOUND IE 14 ENe TSeBGUND2-T 25 
T$4eGQ9XTOVCR40/CQe CL" TOVERWOXTS/ISPOD)XSINCFEEDS 
ICM е GOXxTUVERWO/QxCCS(C(FEEO2/V5 
M] E КЫНАТ (3ЭЗКІМПІ + TCAXLP (4,3KINDJ)> 
Md € «*TsixLP(11,39xXIN0J] + TCAXLP[(12»BKINDJI5 
د‎ VPSHXCPULSTERTND JOY STOMXEPUYS,BK TNIDS 
cUMMENIGUMPUIL THE New VALUE CF THE CONTROL; 
IIT E OOP N OTE OI) +MIXERRLL) +M2xERR (2) 
* M3x«ERR(3)J 
Ir T=zEwTIERCT2/HH+,0000001)+1 AND STAGES=2 THEN 
CUMPUTECHASTeI 
ELSE COMPUTECUASTeOJ 
ПР 150 AND 130 THEN GO TO L3 ENDS 
(4: IE ИМАТИ НЕНА УНИ е ТЕЗРЕИЛХӘМОЕХ ЕНЕМ І%1 ELSE 15 
EEC GCG АРИЈЕ С рт РОСА ТОЈ) = 0) 
R e А[1] + ЧЕ; 
S e-spuCx[431)» С € COSCX[4J)> У e X(305 
COMMENT HERE ARE THE SYSTEM DIFFERENTIAL EQUATIONS: 
UOT E VSS 
Ис ОЙН: 
с е (ше К ООЛО оуу - KXS/R*2j 
ENDO (ХСС + gxSEbNCOFEE)Z VV Ft235 
ӘШИЛЕШІ INEMHNMPSSSECTION THE NEW COAST DURATION IS COMPUTED} 


IF CUMPUTECUAST=1 UR COASTCOMPUTE=1 THEN GO TO LS 


hy 





Иос 2077026; 
t COASICOMPUTLCECUASTCOMPUTE +33 
А(1260; 
СООКО PTS 
Ms le=myvxt7Re2 + FL2) = F(A} 
BKINDEMAXINJEX@ENTIERCT2/HH+ 0090001 D3 
MM1€A(2]xLP( 3»,23KI1N0] « A(3J]XxLP( 4» BKINO JJ 
MM2€A[2]«LPL11»3KXIN0) * A(31xLPf 12» BK IND)J 
MaEBCcAmDPOxLPLI5^»3XINO) * A[3]XxLPÉ15»BKIND13j 
D[1»1] e M41*25 U[(2»2] «e 4M2*25 0[35»3) e MM3x*25 
D(1»2) e Dt2»1J) € MM1xMM2j | 
200153) e ОСЗ ИТТИ ЗУ 
Dt2»3] «€ 0(3»2] * MM2xMN3j 
IF CUFBBITELSIASTEIETTH END ти ба 
COAST € (OLDOCUASTeCUASTO*(QCCMMIXERR( 1)*4 MV2xERRC2] 
*tM4M3XE3XEL31)/LAMDA)5 
COMENT CUAST IS SET TU THE "NEAREST INTEGER DIVISIELE BY 
THE STORAGE INTERVAL» THIS AIDS IN KEEPING INDICES STRAIGHT; 
Kal eee LT ACUAST|T = ТИР О ОРО) 
COAS Wem PESCCCOAST+.907 2) x2; 
ТЕ CURSHeO UR STAGES SI TEN COATES 
WRITES COAST =775,F20, 105 СПА5Т); 


Ló: END SLOPE, 
COMMENT. INTEGRATION PROCEDURE GOES HERE; 


COMMENT MATN PROGRAM BEGINS HERE $ 
COMMENT INITIALIZE BELLS» FLAGS» AND CONSTANTS} 
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ІТЕКе0; Кү Две; FESTE LUNT E00, QUPITeOJ 
RE ée 2,0337; audeas OL XK «€ GOxRE«27 


COMMENT SAMPLETIME IS THE FEEDBACK SAMPLING INTERVAL AND 
HH IS THE DATA STORAGE INTERVAL» HH IS ALSO THE INTEGRATION 
STE PSI AES 
ЗАВЕТ ИЕ 202; Ене 2; 
L1: READCXPL3»0]» XPL4»01],5 TASFF 1»FF2»oWR1» XR2»W1»"2» COAST» TB2» 
LAMDA»STASESI[L4]5 
CONTENT THE INPUT DATA i$ 


xPt3,0) = INITIAL VELOCITY 


i 


XD) s ANN CA. ANGLE 
Т1 = FIRST STAGE BURN TIMES 
COMMENT Реј = FIRS] STAGE THRUST 
FF2 = SECUND STASE THRUST 
aR} = FIRST STAGE FUEL FLOW RATES 
COMMENT wR2 = SECOND STAGE FUEL FLOW RATE 
КЇ LAUNCH WEIGHT 
М2 = SECUND STAGE WEIGHT AFTER SEPARATIONS 
COMMENT COAST = INITIAL CHOICE CF THE COAST DURATION 
132= SECOND STAGE BURN TIME 
LAMDA = COAST WEIGHTING FACTORS 
COMMENT STAGES = NUMBER OF STAGES; 
OLDCOASIECODASTS; 
T2€T1 1 *CULAST; 
TFei2+1c2, 
OLDMAXINDEXeMAXINDEX € ENTIERCTE/HH+.51)5 
ISP1eFF1/WNR1? ISPLEFF2/MR25 


A (2,01 € 0; ХР[Л,0] ¢ XP[4»01/57,29577951315 
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COMMENT 
FOR 


COMMENT 
FOR 
FOR 


COMMENT 


GENERATE OR KEAD IN TRE NOMINAL CONTROL HISTORY; 
PA ЕЕ ПИЕ ИАА ТОЕ DO PHITI)] « 05 


INITIAL CONDITIONS FOR BACKWARD INTEGRATIONS 
je2 STEP 1 UNTIL 22 NU LPLI»0]¢€0;3 
Leslee oo ONT Teor peer Tir ey eT 


ШЕПТЕ Tae NOMINAL TRAJECTORY: 


ROS Urea Os ТЕ» ПИТЕР,» 0» 0» ХР» FUNCT )3 


COMMENT 


PERFORM THE BACKWARD INTEGRATIONS 


L23 ADAFS(2Z,KHH, 0» TF» HH» TFs Or Os LP > LINBAK)J 


COMMENT 


STURE THE ULD VALUES CF THE STATES BEFCRE COMPUTING 


A NEW FORWARD TRAJECTORY, 


FOR 1660 ӨЛЕР 1 UNTIL MAXINDEX 00 LOOPL OLDXPCL»I} e XPILs1D3 


COMMENT 


CEES TIRKE AND LS TBE ALTITUDE CHANGE DUE TO 


PES INES TS "СОИРОТЕО» 


PREHF € ХРГ іх МАХІМОЕХ) 57 | PREVF € XP[ 3»MAXINDEX)J5 


ОБ ИО РЕ МСЕХ 13 


Sheet pd 20, MAXINDEX 15 
Ji3€e€LP[ 22; MAXINOEX]5 

Pe (ieee el PL ios MAX INDEX), 

[КРО Je THVP(L 1» 2] e LPE21» MAXINDEX)5 
ThP(2,2J)eL PL19 »MAXINDEXJS 


INVERTCTMP,251)5 


II EIN HRITEC<"THE CAMEDA TT NATRIX IS SINGULAR"? )5 


GU TU L^ ENDS 

ШОК Ге) ех ыРГ]»1]+,) 1 3® ТИР 2,11; 
ШР НЫ )е)12Х]МР[1,2)+у)13хХТЫР({2›2]; 
VU ¢ SGRTC(KZCPREHF+RED)S 


PREDHFOeTrUPL I2 11X CVU-PREVF Oe TMPC1,2xPREGAMFS 
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PREHFÜEPREHF «*FREDFFOS 
ES O ESO IE оо асо Р КЕНЕС? 
SAT? ГИДО S0 THEEEGIN 
DHF eFREOHF OS 
OVE eSGRTCK/CPREHF * DHF *RE))7FREVES 
0 бАМҒ«-РКЕСАМҒ) ЕМІ)) 

HOUND é BOUNU2 « CUASTCOMPUTE e COMPUTECOAST € 05 
COMMENT COMPUTE THE NEW CONTKXCL. WHILE INTEGRATING THE 
SYSTEM EQUATIONS FORWARD; 

ADAKS( 4s ths Os TF» HH» TF» QO» O» XP» SLOPE) 
COMMENT Tre FUOLLOwING LOGIC CAUSES THE CONTROL FISTCRY 
ШОГО РЕГИ РЕ OUT EVERY EFUOURTN ITERATIONS 

IF ENTIERCITER/42X4 = ITER THEN BEGIN 

IRITECFOR И Т е5 БЕР 5 UNTIIDEEMAXINDLX DO PET 
ЕНІГШЕСІРДСЕЙО ЕМ; 

ТЕ ЕНЕН РО 

OLONAAKiNDEXeMAXINUÜEAX; 

TPCA S ISSUED OOAST ^ THEN 


BEGIN 


COMMENT ADJUST STORAGE LOCATIONS OF PHI AND FLY NEW 
NOMINAL TO PROVIDE PROPER INPUT OF STATES FOR BACKWARD 


INTEGRATION?» 


12611+С04512 
TFel24T32; 
UFDHARTEBEXENAX INDEX 


MAXINDEXeLNTIERCITF/HH* 45131), 
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CENA US А УНА. 0090001 у2 
РОА ЈЕНИ РАИ ин КОТ ЇЇ ноја о ин. Ообосрожу STEP 1 


UNTIL OLOMAXINDEX DO BEGIN 

PE ETL TEPA 

Ie CENTIER UTZ 7HR ОО ОООО =1 ІНРМЕРНІСІЗІР1 60; 
ENDS | 
])e£NIL IER(CT1/HH* 0000001); 
ПИ are AP ELs 15 
Ape cde Hs lie ire nits Tes Os Os Yer FUNCT) > 
FOR rel STEP UNTIL HAZ TNDE F00 
LIPPI AAPP TIS 
END; 

(ОКК ЕК PRINT DUT THE RESULTS ЈЕ TRE ITERATIONS 
NRITECF2); МАТЕ СЕЗ, ОРГ» OVE: DGAMEXx 57.2257 И НИЗ 9909 
WRITECF4»55 | NRITECF3» ХРГІ»МАХІЗГЕКІ -* PREHF» 

ХРГЕЯТФКРЯОСА) + PREVF+ (PREG) MAX INDRA] —PRE G'AMF ) 
x5 FEST 3:95 
WRITE (r5» XP[1l,MAXINDEXJsS XPE3^F AXINDEXJ» 
APLI vASNDEXgXx5!.ec997 789509805 SORTCK/CRE + 
XPCI1» MAXINDEXx 22); 
NRITECFO» COAST), 
COMMENT CUMPUTE THE CORRECTED TERMINAL ALTITUDES 
VJeSQQRT(ÁKA/ CXPLI1o2 MAXINDEXJI«* RED); 
Dem СҮЙ ХРЕЗУМААТИрЕХТ) " 
TMPL1»2]XXPLA» MAXINDEX]J 
ICRA AAA АЈ + DHFOS . 
AUN? ROO. TR CURRECTED TERMINAL ALTITUDES 


ЕШТЕСРТ>НГО)» 
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БАС И t=] THEN GU TO 6) 
COMMENT ГЕСЕ ОКАН МАК ERRORS WITHIN TOLERANCES 
IF Avot PVret VUTAR TS. 14A x VNDEX)))>50 
OR ABSCDGAWFECAPDL4SsMAXINDEXI)) > „02 
Ü4 SIGNCDHFOZSIGNCXPUT^MAXIRNDEXI7"PREHEFO THEN 
BEGIN 
FALL emali+i> 
CUUNTeCOUNT*15 

COMMENT TEST FOR I&/PROVEMENT IN SATISFYING THE TERMINAL 
CONSTRAINTS OR IMPROVEMENT IN THE CORRECTED TERMINAL ALTITUDES 

BP (ABS(OVF )SABSCSSRT(X/CPREHF4RE)) -PREVFO AND 

ABS COGIDO SN O IMPRESA) y" OR (CHED > RREHFO) THEN 

BEGEN 
РРА САОН ЕС EDEN TSV EINEN O ELSE DHEUDS 
530 fU LS END ELSE BEGIN FLAGe0; 
FOX Jeg STEP 41 UNTIL CMAXINDEXeOLDYAXINDEX) DO 
BEGIN COMMENT IF The ITERATION WAS UNSUCCESSFUL 
DISC TD IME URTA AND ATTEMPT ONLE TO SATISFY 
TERMINAL CONSTRAINTS3 

BREIT JEeOLDPHIELTI 

BEINE EULDXPLL, CI 

ENUJ 

CDASTe0LOCUASTS 

T2€Ti+CUAST> 

lh Sizes 

IF JTIER=1 (HEN BEGIN 

MRT e CL PAGE I )3 


CONMENT IF THE FIRST ITERATION IS UNSUCCESSFUL», THE 
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FOLLOWING MESSAGE 15 PRINTED OUT; 
WRITECS"THE TERMINAL CONSTRAINTS HAVE BEEN VIO", 
"LATED EXCESSIVELY BECAUSE" »//»"THE М0", 
"MINAL CONTROL OR THE INITIAL CONDITIO", 
"NS ARE",//."NÜT CLOSE ENOUGH TO THE ", 
"ÜPTIMUM...., GUESS AGAIN">)3G0 Tn L4 
END; 
GU EB ED 
ENU 
ELSE BEGIN FLAGel; DAFe100000/(24FAIL)} COUNT«0J 
COMMENT THE FOLOAINU STATEMENT CONTROLS PROGRAK TERMINATIONS 


Tt DHF S1009 THEN RESIN DHF eDHFOS QUI Te1s GO TO LS 


ENDS 
ENDS 
I» DVFeSQRTCK/CAPLi,»MAXINDEX3*DHF*RE2)-XPC3^,MAXINDEX)5 
ОбАМҒезХРГа»МАХІМДЕХІ?) 
Ol DE 3 
COMMENT .Fl'"AL PRECISION RUN, 
L61 AYAMS( 4» 1» 9» ТЕ, 20» 20, 874, G-45» XP, FUNCT) 5 
CUN MENT ВИА SENSITIVITIES GESIEANÜORE TO ERRORS IN 
TERMINAL CUNSTRAINTS; 
Wiehe cme tere (HP Lise) )> 
COMMENT PUNCH THE JUPTIMUM CONTROL HISTORY ON CARDS; 
hucCRINBECOARDSOAHISTURY, FOR-Leo STEP 1 UNTIL MAXINDEX OO 
ОЧИ 


SMOD; L4s ENO, 
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APPENDIX B. SECOND VARIATION COMPUTER PROGRAM 





BEGIN COMMENT ERNEST C. LUDERS BOX 219 SECOND VARIATION; 
REAL RE»G95,X5, ISP,TF, TOVERNO»NU1»NU2, HH, VO» TE1^ TE2, 
R»V»,TSING»OLODETM^; DETMSOLDTE1»OCLDTE2»HUK»SO T3, T4, 
0105171,» OLDNU2, OLOT) 
REAL SUM» FEE; MASRAT^ $5, C,Q,C1^» C2, C3, CO» INT; 
REAL ARRAY ТЕМРЗ,ТЕМР6ІГ1:31» ТЕМРА» ТЕМР5(1:3,1:231» 
DELX,TEMP1» TEMP2»8TL1:3J3» ^INV»M»NE1?:3»1*:3)» 
COLO а о Ио Јо 
INTEGER | IND» BKIND; 
INTEGER JTERsMAXINSEXs Le ONCE КУ» Ј» STS 
INTEGER BOUNDS 
REAL ARRAY OLDXP»XP,YPLO:4,0:125),PHI>HU,HUVEL0:125), 
FJ RUXLI: 370:1T23J7 ZZ[072403 0:1259 
ELS SLO OS AE ӘН»; 
БАР ые ур) 
БАШ А ТАДА ОУН А ССЦ 18. 12); 
FORMAT FICTON THIS ITERATIONS SKIS SN] zo SEDO 
//»X5» 
"NU? = ОЕ ОО или лахбе TEI SE ЕСО Шр 
KISS TE ЕЕ Ауу 


CEO O AUK “ste COs 11)» 


DEFINE LonpI POR@ let poy sewers LOOP = FOR 261,273 0007, 


puppi S FORS P2 D) 7; 


CUMMENTOSMAERIX INVERSIOMSPROCEOSURECGOES. HERE 5 
PROCEDURE NOMINAL (T>X>DX)5 


VALUE Ts | REAL T» REAL ARRAY X»DX[1)5 


за 





BEGIN 

REAL S2C,MASRAT»FEE»GJ 

INTEGER 135 

Poe NTT TARD J E TIF ICMAXINDEX THEN I+1 ELSE lS 
EE OG R(TOS m. C T 7 titimsoreoPRHI Je); 
REX[1I+RES 

MASRATe1-TOVERWOXT/ISP3 

БСО CERCOS OLAT УСХГ ЗА 

DXLile VxSj 

DXC2)¢ VxC/R> 
DXL33 * (с © GOXTOVERWO/MASRAT)xCOSCFEE) 7 KxS/R*2j 
DX[4] € »K/VXC/F*z + GXSINCFEE)IZ/V «Dx(215 


END МОМІМА(; 


PROCEDURE ВАСК(Т»7,07); 


VALUE Ts REAL Ts REAL ARRAY Z»02[1)5 


BEGIN 

REAL Т8»55,А»|АМ1»| АМ2»| АМЗ,КС»К2,5853,У2»5,С»5Ғ»СҒ, 
МА5ЧАТ»У,1ҺКТ; 

REAL ARRAY ҒХ,НХХЭ»Ғ,0,55,М»М»5МЭЕТМ»ҒМ»0М(1:3,1:31,СС, 
DsMTISNTCeHXULCL? 313 

INTEGER Ej 

LABEL EN; 

TBeTF-T;  I*ENTIEPCTB/HH*,000000001); 

IF ISO THEN I є 0; MASRATe1-TOVERWHOXTB/ISPj 


J © IF J«MAXINDEX THEN I+1 ELSE I; INT *TB/HH ~= l} 


WN 
WN 





R € (UexPL1» 13) * INTxCXPL1»J] = U) + RE; 
ccu SAMO ENIXCXPIE3) J] = Шр; 

SEES TIO CE TNTCCXPES, JJ 999 
C ceu се) ) t INTXCXPLA&4»J1 = W))> 
SFe SINC CUERHITCII) + INTXC(PHIDLJS - 0/92; 


Свен ос ноше Ре тј) + INIXCPHICL JJ]  0U))> 


COMMENT COMPUTE PARTIAL DERIVATIVES OF SYSTEM EQUATIONS 


WITE RESPECT 10, STA ЊОЈ 


РХ ТА аса ЕСИР ИЕГИ Зеус; 


ЕХСГ2,11«2хКх5/Ққ»3; ҒҰХ(2,2160; Fx[l2,3]Je"“KxC/R*253 


FXL3,1Je-VxC/R*2«2xKxC/CVxRx*3);j 
FXE3»2)eC/R*KxC/(VXR2)*2*7* GOXTOVERANOXSF/CV*2xMV ASRATOS 


FX[3)3]e-VxXS/R+KxS/(VxXR*2)5 


COMMENT COMPUTE PARTIAL DEFIVATIVES OF SYSTEM EQUATIONS 


NITH RESPECT T9 CONTROL; 
FUL1»11«05 
FUE2» Ile-CA*GOXTOVER^O/MASRAT)xSF;j 
FUE3, I)eA/VXCF 5 


COMMENT COMPUTE DERIVATIVES OF THE HAMILTON] AN? 


LAM1eZL[22)5 LAM2eZL23)7 LAM3eZ[24); 
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BUI TJeLAM2XFUL2> LI+LAM3XFUL3HI 
HUULIJe-LAM2xAxCF-LA"V3xAxSF/V = HUK5 


COMMENT CHUCK HUU 2 03 


IF HJULT320 THEN 
IF TSING=0 THEN 
BEGIN 
TSINGETR+HHx5; 
WATTECC"HUU 2 O AT TIME s "»F204.10»5» TB2J 
YRITEC<"HUU = “,£20,10>»HUUCITIS 
If TSING2TF THEN 
REGIN 
TSING*0j 
TSINGeti/TSING 
ЕМО? 
GO TO EN 
ENO 
ELSE GO TO EN; 
HUXC1»IJeHUXC3» 1Jeuxut131eHXxUC3)*0; 
HXUL21eHUX(2» IJ "LAM3xFUL3^ I]/V5$ 
KSeKxS; ксекхС» Р2«<5хК; Угеухуз БЗекгехр; 
HXXL1»11€e€-6XLAM2ZxKS/R2*2*LAM3X( *6xKC/(R2xR2xV)*2xVxC/R3)j 
НЕ Хе ен ха Асо ХК МохХЕЗ) + С/82); 
HXX(1»3J€«nXXL3^»1)€*€2xLAV2xK C/R3-4LAM3x(-"2xKS/(R3xV) *VxS/R2)j 
HXX(2»323J]e-LAM3x(-2xAxSF/(V2xV) -2XK C/ CR2xV2xV))j 
HXX(2»3JeHXX( 3» 2) *LAM|1xC + LAM3XC=KS/CR2XV2)=S/R)} 


HXX[3»3Je*LAM1XVxS * LAM2xKS/R2 * LAM3XxC(KC/CR2xV)-VXC/R)j 
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Kxe05 
LOOPL LOIPJ 
BEGIN 
F(L» JJeFXEL» JJ7FUCL^ IJ3XHUXCJ2 I37HUUE I 27 
ОГ ll Js 1J3/78UUCTE3) 
SStL»JJeHXXEL^ JJ-HXUECLJ. xHXUCJ] /HUUCIJ; 
MELAJJeZECKK*KK«1235. NEL» JJ2eZEKK«91; 
COMMENT STORE F AND SS MATRICES FOR USE IN PROCEDURE NEWLAMDAS 
IF I2100 THEN 
BEGIN 
FLCL^J^ IJeF tto 225 
SLEL^,J»^ 1JeSStL 5» J2; 
END; 
END; 
IF I<MAXINDEX THEN 
DETM-ME1,1Jx(M[2,23xM[ 3^» 3)»M(2,3)xME 3,2)])* 
М(1>21Х(М(2,31хХМ3»11М(2,11хМ(3,319% 


М(1»31х(м"72,11хмМ(3>21"м(2,21хМ(3»1109; 


COMMENT CHECK FOR CHANGE IN SIGN ОР DETERMINANT ОР М} 


IF I«CMAXINDEX-1) AND 
SIGNCOLODETMO)ZSIGNCDETHM) THEN 
IF TSINGz0 THEN BEGIN 
TSINGETB+HHXS5 
ЕКІТЕС<"ТНЕ DETERMINANT OF M CHANGED SIGN AT T="» 
220.10>»78)2 60 ТО ЕМ 


ELPEDGUNSE GO TO EN: 
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AeHULII/HUULI)Z 
POOP GEGEN 
CCECJJe-FULJ^ IJx45 
DLJJeHXULJ) x As 
END; 
OLCDETMeOETM; 
00Р ВОВЕ ЕС 
ЕмМЕЈ,1 7607 
GNLJ,L3*05 
SMLJ,L]JFO3 
ЕТМЕЈ, 1 60) 
FOR «K*1»2,3 DO REGIN 
СОЕ ЕЕЕ ЦЕ ОКК 1) 
QNCJ»L Jeg NC Jo L 1e GEJ^KK2x МЕККУ 77 
FTNEJ,UJ*F TNCJLJ*F UEK^ JJXNETKK,L35 
МЕЧА ез Ол ]+SSCJsKK]x< MIKKsL}S 
END; END» 
KK eos (UCP EI OLE F NIC) «O72 bte1 eo, 
COURE (Uses “SEG 
DZ CKKEKK 41) Je МОША ЈОВО >» Ј Ј> 
ПАКУЈЕ УМОВА ЕРЕ Уа Јо 
ВИСЕ е1 +2107 ОБЕ ТЕЕ а 24 1 
MTOLLJIEMTDELI] + МЕЈ, 1х0ЕЈ7; 
NIESELEIeNTCLLI + NEJALIXCOLIYS 
END; 
LOOPL РАДЕ Ва ИЕЛИ МТСГ Зх 
EN’ END BACK; 


БАРЕВО ВА СОМЕ 1 rx eX) 


2 





VALUE Ts REAL Ts REAL ARRAY X»DxX[1J5 
BEGIN 


РАВЕС пи Ре 2 3» 24 


IF T20NCEXx,99999699 THEN 
BEGIN 
IF BOUND*O THEN 
BEGIN 
INDEENTIERCONCE/HHA+,000000135 
о ош ко 2 
END; 
L1 INDeENTIERCONCE/HH*,0000001)*15; 
ОМСЕ«ОЧСЕЗНН; 
(2! BXINNEMAXINDEX=1ND5 IF IND2MAXINDEX = 7 THEN GO TO L4; 
KXxe05 
LOOPT 3EGIN 
ВТО tier ak inva, 
LOOPJ BEGIN 
MINVCI^JJeMCI»JJeZZUCCKK KK t i2)» BKIND]J 
NECI» J3*22U0KK*9» BKIND) 
END»; 
END; 
INVERTCMINV,? 3,125; 
HEMOS THEN "RITEC«"M IS SINGULAR AT T = "F6,2»»155 


COMMENT CO VTE INITIAL CONDITIONS TOR PROCEDURE NEWLAMOAS 


IF 1ND=100 THEN 
LOOPI SEGIN 


LOOPJ BEGIN 
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TEMP4OCT,JIJ©MINVOI,JI3 
eee ails Pe NCI,J)3 
END» 
TEMP6LIJ*BIIJ; 
END; 
IF IND=102 THEN 
BEGIN 
DESSEN Ne) OLNXPL1>,1007>, 
DELX[2]eXP[ 3; 100]F0LDXP[3»100]; 
Dele Seno foe OCl=-CUpXPl 4,100); 
ENORI BEGIN 
SUMeO; 
EDOPJ SUMESUMYIEMPS[LJ,IJxDEUXUUJIS 
TEMPILIJeSUM+TEMP6LIIS 
END> 
LOOPY BEGIN 
SUMeOs 
LOOPJ SUMeSUM«* TEMPAC J,IJXTEMP1UCJJ; 
2РГІ»01«50/М; 
END; 
END; 
COMMENT COMPUTE COEFFICIENTS FOR COMPUTING NEW CONTROL; 
L9GPI REGIN 
SUMeO5 
LOOPY SUM*SUMN4FUC JU, INOIXMINVE Is J)5 
TEMPILIJeSUM 
END; 


LOOPI BEGIN 


29 





SUMeO0; 


, 


mOUSMESUMCSUMTTEMPILJIXNET» JJ)? 
ТЕМРӘГІ7650/М; 
END; 
(42 IF InD >MAXINDEX THEN GO ТО L33 
TEMP3L1)Je0LOXPL 1, IND13 
TEMP3[2]eOLDXP[3,INDJ5 
TEMP3E 3) e9LDXP[ 4 IND); 
13eT4e0;5 
LOOPJ BEGIN 
COC EDI SOHO - TEMPOPULJI)/HUUELIND)S 
T3*-T3+TEMPILJIXZZLIJ+18,BKINDJ3 
END; 
T3*CT3+KULINDII/HUULINDI 
BOBBY  TqgeTrTaecCOCUTND)x TENES JO 
Sm IND )EPRICINDI= 749 T3S 
NRITEC«"INDEX = ",14>,1 0); 
WR ME pps" ESOO) CEE ЛЕССИ О "СЭ = ",E20.105 
ИИ Боб, О POR IES TES UNTIL. 4-DO cCcCLIS INDI, 
ИКТТЕ(ТЗ,Т4)) 
IF BOUND=0 THEN 
BEGIN 
BUUNNeRBOUND+15 
uS mer 
END; 
END; 
L3: I © ENTIERCĪ/HH); Je IF I<MAXINDEX THEN I+1 ELSE 15 


ТЕ ТАБЫТ ТНЕМ 
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BEGIN 
OLOTeT; 
INTeT/HH71; 
616(0«СПСГТІ,ГТӘФІМТХІСПІ1,,/1-0); 
C2e(CQecOF 2^» I139)« INTxCCOL2» 21-02; 
C3et(Q0e20L3,122*INTxCCOL3, 12-02; 
C4etcQaecaL 4,12) * INTxCCOL 4, 1-025; 

END; 

КЕКерРНТ(Т/НН]#С1хХ[1]+С62хХЇ{3)+С63хХх{[47]+С4; 

ReXLi)+REF MASRATe1-TOVERNOXT/I1ISPS5 

+5 INSIDE BSD VeX£3)3 

OX ts US, 

DXt21e VxC/RJ 

Dxl3Je ( 9eGOQxTOVERWO/MASRATOxCOSC(FEED) 7 KxS/R*27 

DX(4)e7kK/VxC/R*2 «* GxSINCFEED2/V + DX[2)5 


END CONTROL; 


COMMENT | INTEGRATION PROCEDURE GOES HERE; 
PROCEDURE NEWLAMDA(T»LAM^DLAM)J 
МЛАДЕ РА Ts REAL ARRAY LAM»DLAM[1JI5 
BEGIN | 
REAL д» ІМТ»0; 
REAL ARRAY D»DELX»^SDX ,FTDLEL1:32,5DELXC1:3,100:110) 
МЕОИ 
INTEGER IND»KKsJJ5 
BES kl lL 
IF TZQONCEx,990090999 THEN 


BEGIN 
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IF BOUND=O THEN 
BEGIN 
INSeENTIERCONCE /HH*,0000001 55; 
GO TOML2 
ЕМО? 
(212 INCC€ENTIERCONCE/HH«,0000001)2*15 
ONCE e-ONCE * HH; tee i Nb ow AY INDEX THEN GO ТО LU 
82: DELXLIIEXPLI>INDITOLDOXPL1- INDIS 
ПЕЛРХАРОН МРИРО ОФ -ЦШКОХРГЗ»ТАЫрТ» 
DELXC3)exPEC 4A, IND)"OLDXPL 4» INDJ25 
Ae HUT INDI/HUU[INDI5 
LOO?I BEGIN 
DEIIeHUXLI,INDIXA; 
SDELX( 1, INDJ+OS 
КООРУ SDECXCISINDI]eSNELXLI»INDI+SL[IsJsINDIX 
DELXLJD 
SDELX[I; TNO SDEEXET , 1 NDS > 
END? 
IF BOUND=0 THEN 
REGIN 
BOUND BOUND +13 
GO TO L1 
END» 
END; 
MER KXeENTIERCT/HRIS JJeTF KK<MAXINDEX THEN KK+1 ELSE KK} 
| КЕ ТУПАК» 
LOOPI BEGIN 


Sie Ol Сес ASAS IELNTZESDELXCTI>. JJI 0); 
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Ca ASE COEL TJ) KK II+INTXCEL[l+ Jr 1110)5 
END; | 
LOOPI BEGIN 

FTOLCIJe«0; 

Pies FUSE LTIIeRIDELCIJ+FLIE YU, 13 XLAMC Jl? 

О АМС12е -ӘБЕІТТЕТПІСІ?» 
END; 


END NEWLAMDAS 


COMMENT MAJN PROGRAM BEGINS HERE; 
RECe2,00€7; GO*32.173 KeGOxRE+2;  ITEReO; ISP+3005 
0.017.900; 
КЕАОСХРЕ З, о] XPE 4, 03)» TF, TOVERWO^»;NU1» NU2» HH, HUK 5 
ХРГ1»61еХР(2»01е0; 
ХРГЛУ ФИ ЕХРГДЫОТ/Ф 22957795131; 
BAKDNBEXeENTIERCTFZHH*45127 
READ@MESTORY»,FOR JeO STEP 1 UNTIL MAXINDEX DO PHILIJIS 
ARITECHISTORY»FOR 10 STEP 1 UNTIL MAXINDEX DO FHILIIIZ 
COMMENT INTEGRATE SYSTEM EQUATIONS WITH NOMINAL CONTROL; 
ADAMSCAsHH, Os TF» HH» ТҒ» 0, 0») XP» NOMINAL); 
VOeSQGRTCK/CXPCI1»MAXINDEXI *RE2 25 
Nu2e*€NU?2/(2xV0)» 
ПЕТЕ КРЕЧ MAXI 2 
TE2€-VO . -XPL3»^MAXINDEX] 3 
BIG ІР ТРЕБАО ANO AEE ODED E GD EDED ASSCOLDTE2)«48SCTE25 
_ АМО OLOXPC1»VMVAXINDEX J»2XPLI1» MAXINDEX)J THEN 
BEGIN 


ТЕ1ТЕ1/2; TE2e 72272 HUKeHUKx105 
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END ELSE HUKeHUK/23 
COMMENT COMPUTE INITIAL CONDITIONS FSR BACKWARD INTEGRATIONS 


ReXP[1»MAXINDEXIS+RES VexP[3>MAXINDEX]I5 

БІН ТӘТЕ UNTIL 18 DO ZZ O EOS 

ZIEH Dear VE 

22[4»0) «€ K/R*2) 

2711001 «€ -üxVxKxNU2/R*3j 

217{12,0] e -77(4,0]; 

22113,0] е -2xKXNU2/R*25 

LZ2U15»0) e«22t1»035 

27{17»0) ée 1; 

ASIS OT > PELO 

22(20»016 "TE, 

222170) £03 

IF ITER=SO THEN BEGIN 

22022»01e1«KxNU2/R*2j 

27(23»01в6 2xVxNU2j 

22l 2O je NUL, 

ED ELSE BEGIN 
BOUNDeO; 
ONCE «200; 
ADAMSC3,HH,200»5 TF, TF» TF» 87/1, 8*5, ZP, NEWLAMDA); 
LOUP TOME  Ү+772 [22,0% 
горе оледи тео зх Од 
ZzB2u, Veen) 7202005 

END; 

Бәкей» - ONCES<OS 


DETMEOS 
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COMMENT | INTEGRATE M,N»B, AND LAMBDA EQUATIONS BACKWARD; 
ADAMS 24> Whe Os Ths НЫН» 10» 0 >0 » ZZ»BACK);j 
COMMENT STORE OLD VALUES DF STATES BEFORE RUNNING NEW TRAJECTORY; 
FOR I«0 STEP 1 UNTIL MAXINDEX DO 
РОК belua, 3,5 DENIED eXP Eb» T5 
COMMENT IF DET(M) CHANSED SIGN THEN A CONJUGATE POINT EXISTS 
OR IF HUU WENT POSITIVE DURING THE INTEGRATION THE LEGENDRE 
CONDITIONIS NOT SATISFIED, ONG IS SET TO THE TIME AT 
WHICH EITHER OCCURRED AND THE FORWARD INTEGRATION IS STARTED 


AT THIS POINT RATHER THAN ZEROS 


ОМСЕв0; 
BOUND«O, 
IF TSING=N THEN ADAMSCasHR»0s TFsHH» TF>0s0» XP, CENTROL) 
ELSE BEGIN 
IF TSING>TF THEN BEGIN 
WRITEC<"HUU WENT POSITIVE OR A CONJUGA"» 
“ТЕ РОТ ОТ BCC ESD IOO CLOSE TC ITHE 09 
"END OF THE T3AJECTORY'"»):60 TC L2 END; 
ONCE*ENTIERCTSING*.1)2LeENTIERCTSING/HH* ,1)7 
FOR 26123,4 р0 ҮРІ./,01«ХхРГ2,12) 
АрАМ5(4»НН»ТӘТАС»ТЕ»НП»НН,2»0>ҮР»СПМТЕПІ 0) 
FUR I«L STEP 1 UNTIL MAXINDEX 00 
FOR EE CPI IeYP CJ» T« 3» 
ENDS 
НЕТТЕСЕОВ 160 STEP 1 UNTIL MAXINDEX DO PHILI))3 


\0©#5О0КТ(‹(К/(ХР{[1,МАХ]МОЕА]+КЕ)); 
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OLOTE1leTE1; EEDTEZe ies 
TE1e7XP[LA4, MAXINDEX)1; 
ПЕ Сет Им КХР ИЗ НА ИЕ 7777; 
COMMENT IF IMPROVEMENT IN TERMINAL ERRORS AND PAYOFF IS 
LESS VEN ET STLON CUSERS CHOICE) THEN STOP THE PROGRAMS; 
lf eee Г ОТЕТ ЕЛ АЗ-СВІПрТЕ2-1Е2)<1 
AND ABSCOLDXP[1>MAXINDEXI=XPT1»MAXINDEX])S500 THEN 
BEGIN 
WRITECPOR РО о БЕО 2 UNTIL MAXTNDEX DO PHITT I); 
GO TO L2 
END ELSE 
SEGIN 
ITEReITER+15 
WRITECF1,NUl^NU2,TEX,.» TED» HUX 25 
GO TO L1 
END; 


|2, ЕМО, 
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